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The formalism of matrix product states is used to perform a numerical study of (l+l)-dimensional 
QED - also known as the (massive) Schwinger model - in the presence of an external static ‘quark’ 
and ‘antiquark’. We obtain a detailed picture of the transition from the confining state at short 
interquark distances to the broken-string ‘hadronized’ state at large distances and this for a wide 
range of couplings, recovering the predicted behavior in both the weak- and strong-coupling limit 
of the continuum theory. In addition to the relevant local observables like charge and electric field, 
we compute the (bipartite) entanglement entropy and show that subtraction of its vacuum value 
results in a UV-finite quantity. We find that both string formation and string breaking leave a clear 
imprint on the resulting entropy profile. Finally, we also study the case of fractional probe charges, 
simulating for the first time the phenomenon of partial string breaking. 


I. INTRODUCTION 

The confinement of color charge in quantum chromody¬ 
namics (QCD) is one of the beautiful key mechanisms 
of the Standard Model. Focusing on the static aspect 
of confinement, one can probe the theory with a heavy 
quark-antiquark (qq) pair and examine how the modified 
ground state evolves as a function of the interquark dis¬ 
tance [T]. For small distances a color electric flux tube 
forms between the pair, resulting in a static potential 
(i.e. the surplus energy of the modified ground state) 
that grows linearly with the distance. This flux tube 
can therefore be conveniently modeled by an interquark 
string with a certain string tension. One can then de¬ 
scribe a heavy quarkonium state as a qq pair that is kept 
together by this confining string. However, there exists 
a critical distance at which the string breaks. Beyond 
this distance the flux tube disappears and the potential 
flattens out to a constant. At this point it has become 
energetically favorable to excite light particles out of the 
vacuum that completely screen both the probe quark and 
antiquark, leading to two isolated color singlets. In a dy¬ 
namical setting these would then be the two freely prop¬ 
agating jets of hadrons that emerge as the final product 
of some particle collision. 

This phenomenological picture is corroborated both by 
experiment and theoretical work. At the computational 
level, the static potential has been studied extensively 
over the years with lattice QCD. The linearly rising con¬ 
fining interquark potential has been obtained, both in 
the quenched case that excludes dynamical light 

quark degrees of freedom and in the unquenched case 
p^HTS] that includes these degrees of freedom. In the 
latter case, where the dynamical quarks can screen the 
heavy probe charges, the phenomenon of string breaking 
has also been observed m as an asymptotic flattening of 
the calculated potential. Nevertheless our understanding 
of confinement is incomplete: the Euclidean space-time 
lattice Monte Carlo simulations cannot access the real¬ 
time aspects of the dynamical string formation and string 


breaking. Furthermore, even in the static case, it is not 
settled yet if one can fully describe the confinement mech¬ 
anism - specifically the nonperturbative string formation 
- in terms of (semi-) local degrees of freedom (e.g., center 
vortices, magnetic monopoles) [2Qll22] . 

In this paper, we study how confinement and string 
breaking show up in the Hamiltonian setup, as opposed 
to the Euclidean path integral setup of lattice Monte 
Carlo. We do this for the simplest nontrivial quantum 
gauge field theory: (I+I)-dimensional quantum electro¬ 
dynamics (QED 2 ), also known as the Schwinger model 
[23] . The Schwinger model has a long-standing tradi¬ 
tion as toy model for QCD, sharing its confining and 
chiral symmetry-breaking properties. (We will therefore 
often refer to ‘quark’ and ‘antiquark’ both for the exter¬ 
nal probe charges and for the light dynamical fermions.) 
But notice that in the future the significance of QED 2 
could go well beyond being a toy model, as QED 2 or 
QED 2 -like theories might be realized effectively by quan¬ 
tum simulators [24H43] . 

An important difference with QCD is that QED 2 al¬ 
ready exhibits confinement at the perturbative level, as 
the Coulomb potential is linear in (I+I) dimensions. 
Furthermore, the theory can also be solved in a strong¬ 
coupling expansion, via bosonization [44]. We make ex¬ 
tensive use of both the strong- and weak-coupling results 
in the analysis of our numerical results. Our simulations 
of the lattice Hamiltonian are performed close to the con¬ 
tinuum limit, indeed allowing for a quantitative check 
against these analytic continuum results in the appropri¬ 
ate regimes. 

Specifically, we simulate the modified vacuum struc¬ 
ture in the presence of two probe charges and this for 
different distances and values of the charges. As we will 
show, already in this static case the Hamiltonian simula¬ 
tions give a complementary view on the confining prop¬ 
erties of the theory. At the practical level, the direct 
access to the quantum state allows for a relatively easy 
calculation of all local observables. In this way we could 
not only extract the static interquark potential, but also 
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for instance determine the detailed spatial profile of the 
electric string or the precise charge distribution of the 
light fermions around the probe charges. At a more fun¬ 
damental level, our tensor network state simulations (see 
below) allow for a direct calculation of the entanglement 
entropy between different regions. In the past decade it 
has become clear that entanglement entropy is a very use¬ 
ful quantity for the characterization of quantum many- 
body systems and quantum field theories [45], in partic¬ 
ular also for the investigation of the confining properties 
of gauge theories [46|. In this context the entanglement 
entropy is typically calculated either from the dual ge¬ 
ometry in tire AdS/CFT approach [46l ST] through the 
Ryu-Takayanagi conjecture [48|, or from lattice Monte 
Carlo simulations [49ll5T] through the replica trick, al¬ 
lowing for calculation of the discrete Renyi entropies. In 
contrast, tensor network state simulations give access to 
the full Schmidt spectrum of the state. The Schmidt 
spectrum {Aq,} follows from the Schmidt decomposition: 
if |T) G is a state belonging to the tensor prod¬ 

uct of the two Hilbert spaces 1-La and then one can 
write 


i^) = Ev^ 


a=l 


( 1 . 1 ) 


with d < 


max(dim {Ha ), dim{l-L b)), 


e Ha and 


1-Lb orthonormal unit vectors and Aq,, called 


the Schmidt values, non-negative numbers that sum to 
one. From the Schmidt values one can calculate all Renyi 
entropies, including the von Neumann entropy. In our 
simulations we find that subtraction of the vacuum en¬ 
tropy results in a UV-finite entanglement (von Neumann) 
entropy and that both the string formation and string 
breaking leave characteristic imprints on this renormal¬ 
ized entropy. 

As was mentioned in the previous paragraph, we use 
the general formalism of tensor network states (TNS) [52l 
[53] for our simulations. Although the TNS formalism 
has been mainly developed in the context of condensed 
matter physics, it is actually a universal method in the 
same way that the Feynman diagrammatic approach has 
a universal character. The latter applies whenever the 
interactions are weak, whereas the TNS method applies 
whenever the interactions are local. It is in fact precisely 
the entanglement structure of low-energy states for local 
systems, captured by the so-called area law [54|, which 
lies at the root of the TNS description. 

For one spatial dimension, the most widely used 
TNS go by the name of Matrix product states (MPS) 
[551 El]. Recently, different applications of MPS on 
(l+l)-dimensional gauge theories, have demonstrated its 
potential in the context of gauge theories. In [57ll59] the 
MPS formalism was used for the numerical simulation of 
nonequilibrium physics, but static properties [381I571I59F 
and finite temperature properties [68l [69] have also 
been studied. Notably [59] simulated string breaking for 


probe charges in a S'!/(2) quantum link lattice model. In 
higher dimensions, the TNS formalism is at present less 
developed, nevertheless some first promising results have 
appeared for (2-j-l)-dimensional gauge theories [70117^ . 


In the next section we discuss the starting point of our 
simulations, introducing both the relevant lattice Hamil¬ 
tonian in the presence of probe charges and the appro¬ 
priate form of MPS that is dictated by gauge invariance. 
We then first consider the asymptotic cas e of two (frac¬ 
tional) probe charges at infinity in section HI In section 


JVjwe consider finite interquark distances and study how 
the ground state evolves as a function of this distance. 
We distinguish three different cases. First we consider 
the strong-coupling limit. This is a special case, since in 
this limit the interquark string never forms and all probe 
charges, fractional or integer, are screened asymptoti¬ 
cally. We then go away from the strong-coupling limit, 
considering first the case of unit probe charges. In this 
case we clearly observe the transition from a string state 
at short interquark distances to a broken-string two me¬ 
son state at large distances. Then, in addition to unit 
probe charges we also consider fractional probe charges, 
simulating for the first time the phenomenon of partial 
string breaking, with probe charges that get only par¬ 
tially screened. Finally in section |V] we present our con¬ 
clusions. Technical details on our MPS simulations and 
on some perturbative weak-coupling calculations can be 
found in the appendices. 


II. SETUP 


A. Hamiltonian and gauge symmetry 


The Schwinger model is (l+l)-dimensional QED with 
one fermion flavor. We start from the Lagrangian density 
in the continuum: 


£ = ^ + gA^) -m)^- . (2.1) 


One then performs a Hamiltonian quantization in the 
timelike axial gauge (Aq = 0), which can be turned 
into a lattice system by the Kogut-Susskind spatial dis¬ 
cretization m- The two-component fermions are sited 
on a staggered lattice. These fermionic degrees of free¬ 
dom can be converted to spin-1/2 degrees of freedom 
by a Jordan-Wigner transformation with the eigenvec¬ 
tors : Sn G { — 1,1}} of (Jz{n) as basis of the lo¬ 

cal Hilbert space at site n. The compact gauge fields 
0{n) = agAi (n), live on the links between the sites. Their 
conjugate momenta E{n), with [0{n)^ E{n')] = 
correspond to the electric field. The commutation rela¬ 
tion determines the spectrum of E{n) up to a constant: 
E{n)lg = a{n) , with a{n) G M corresponding to the 
background electric field at link n and p G Z. 

In this formulation the gauged spin Hamiltonian de¬ 
rived from the Lagrangian density (2.1) reads (see 
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for more details): 


and we finally find the Hamiltonian m- 


H = 


A (E 3E(n)» + X:(-1)V. W 




2y/x \ ^ g‘^ 

^ \nez ^ 

-\-x + 1) + h.c.)^ (2.2) 

nEZ ' 


where = (1/2) ± icjy) are the ladder operators. 
Here we have introduced the parameter x as the inverse 
lattice spacing in units of x = l/{g‘^o?). The con¬ 
tinuum limit will then correspond to x ^ oo. Note the 
different second (mass) term in the Hamiltonian for even 
and odd sites which originates from the staggered for¬ 
mulation of the fermions. In this formulation the odd 
sites are reserved for the charge —g ‘quarks’, where spin- 
up, s = +1, corresponds to an unoccupied site and spin- 
down, 5 = — 1, corresponds to an occupied site. The even 
sites are reserved for the charge -\-g ‘antiquarks’ where 
now conversely spin-up corresponds to an occupied site 
and spin-down to an occupied site. 

In the time-like axial gauge the Hamiltonian is still in¬ 
variant under the residual time-independent local gauge 
transformations generated by: 

gG{n) =E{n) - E{n - 1) - |(cT,(n) + (-1)"). (2.3) 

As a consequence, if we restrict ourselves to physi¬ 
cal gauge-invariant operators O, with [0,G(n)] = 0, 
the Hilbert space decomposes into dynamically discon¬ 
nected superselection sectors, corresponding to the dif¬ 
ferent eigenvalues of G{n). In the absence of any back¬ 
ground charge, the physical sector then corresponds to 
the G{n) = 0 sector. Imposing this condition (for every 
n) on the physical states is also referred to as the Gauss 
law constraint, as this is indeed the discretized version 
of dzE — p = 0, where p is the charge density of the 
dynamical fermions. 

The other superselection sectors correspond to states 
with background charges. Specifically, if we want to con¬ 
sider two probe charges, one with charge —gQ at site 0 
and one with opposite charge ^gQ at site k, we have to 
restrict ourselves to the sector: 


gG{n) = gQ{6n,o - Sn,k) • (2.4) 

Notice that we will consider both integer and noninteger 
(fractional) charges Q. 

As in the continuum case we can absorb the probe 
charges into a background electric field string that con¬ 
nects the two sites. This amounts to the substitution 
E{n) = g[L{n) +a(n)] where a{n) is only nonzero in be¬ 
tween the sites: a{n) = —(30(0 < n < k); and L{n) has 
an integer spectrum: L{n) = p G Z. In terms of L{n) the 
Gauss constraint now reads: 


H = ^ [ y][E(n) + a{n)f + y](-l)"(T^(n) 

^ \nEZ nEZ 

^(cr+(n)e*^^’^V“(n + I) + h.c.)^ , (2.6) 

nEZ ^ 


in accordance with the continuum result of m- 

In the following sections we will obtain ground-state 
approximations of this Hamiltonian, for different values 
ofmjg^ different values of the probe charge Q and differ¬ 
ent distances L = k/y^ (in physical units p = I) of the 
charge pair, all this for different lattice spacings l|^/x, 
focusing on the continuum limit x ^ oo. 

An important point regarding the continuum limit is 
that the ground-state energy of the Schwinger model is 
UV divergent but that this UV divergence does not de¬ 
pend on the background field a{n ). If we write Sp = 2Neo 
(with N = |Z|) for the ground-state energy of ( |2.6[ ) with 
zero background field a{n) = 0, we have \/xeo —xj'K 
for the energy density in the x ^ oo limit [79]. For 
the modified ground-state energy in the presence of the 
probe charge gQ pair at distance L we can then write 
Sq{L) = Vq{L) + fo, where the potential Vq{L) is now 
UV-finite. Notice that Vq{L) will also be IR {N oc) 
finite (for finite L). 


B. Gauge-invariant MPS 

Consider now the lattice spin-gauge system ( |2.6[ ) on 2N 
sites. On site n the matter fields are represented by the 
spin operators with basis {|sn)^ : Sn G { — 1,1}}. The 
gauge fields live on the links, and on link n their Hilbert 
space is spanned by the eigenkets {\pn)n ' Pn G Z} of the 
angular operator L{n). But notice that for our numerical 
scheme, we retain only a finite range: < Pn ^ 

Pmaxi'^ + !)• We will address the issue of which values 
to take for Pmin{'^ + 1) and Pmax{'^ + 1) later in this 
subsection. Furthermore, it will be convenient to block 
site n and link n into one effective site with local Hilbert 
space spanned by {\sn,Pn)n}' Writing tZri = {sn,Pn) we 
introduce the multi-index 

« = ((si,Pl), (S2,P2), • • • , (S2JV,P2JV))= 

With these notations we have that the effective site n 
is spanned by {\i^n)n}' Therefore the Hilbert space of 
the full system of 2N sites and 2N links, which is the 
tensor product of the local Hilbert spaces, has basis 
{\k) = \k,i)i ... \i^ 2 n) 2 n} ^ general state |T) is thus 

a linear combination of these |k:): 


cTzin) + (-1)" 


(2.5) 


with basis coefficients G C. 


G{n) = L{n) — L{n — I) 


2 





4 


A general MPS |4^(A)) now assumes a specific form for 
the basis coefficients m- 

|^(A)) = ^ 4A., (1)^.,(2). .. {2N)vr |k) , 

K, 

(2.7) 

where (n) is a complex D{n) x D{n-\-l) matrix with 
components [A^^{n)]ap and where vl ^ ,vr G 

^D(2Ar+i)xi boundary vectors. The MPS ansatz thus 
associates with each site n and every local basis state 
\i^n)n = kn,Pn)n ^ matrix A^^{n) = As^^p^{n). The 
indices a and /3 are referred to as virtual indices, and 
D = max^[D(n)] is called the bond dimension. 

To better understand the role of the bond dimension 
in MPS simula tions it is useful to consider the Schmidt 
decomposition (O) with respect to the bipartition of the 
lattice consisting of the two regions Ai{n) = Z[l,... , n] 
and A 2 {n) = Z[n + 1,..., 2A^] [56] : 


D(n+1) 

mA))= 


( 2 . 8 ) 


a=l 


Here (resp. are orthonormal unit vec¬ 

tors living in the tensor product of the local Hilbert 
spaces belonging to the region Ai{n) (resp. A 2 {ti)) and 
AQ,(n), called the Schmidt values, are non-negative num¬ 
bers that sum to one. O ne ca n easily deduce that for a 
general MPS of the form ( |2.7| ) at most D{n^ l) Sc hmidt 
values will be nonzero (for the cut at site n (2.8)). We 
refer to appendix and for the computation of the 
Schmidt values for the specific case of our simulations 
and to [56l [80] for the general case. We thus see that 
taking a finite bond dimension for the MPS corresponds 
to a truncation in the Schmidt spectrum of a state. The 
success of MPS is then explained by the fact that ground 
states of local gapped Hamiltonians can indeed be ap¬ 
proximated very efficiently in D [ST] and that the com¬ 
putation time for expectation values of local observables 
scales only with allowing for reliable simulations on 
an ordinary desktop. 

Another advantage of MPS simulations is that one 
can work directly in the thermodynamic limit N ^ oo 
[801 [82], bypassing any possible finite-size artifacts. In 
the following, we work in this limit. In section [IH| where 
the Hamiltonian is invariant under translations (over 
two sites), the tensors A ^^ (n ) depend only on th e pa rity 
of the site n, see eq. (3.1). While in section [TV] the 


MPS ansatz is not translational invariant in the bulk, 
see eq. ( |4.1| ). In that case the tensors will be fixed 
asymptotically (n ^ 1) to their ground-state value, 
anticipating that we approach the translational invariant 
ground state of the zero-background Hamiltonian. In 
both cases the MPS ansatz depends on a finite number 
of parameters. Finally, we note that, in the thermody¬ 
namic limit, the expectation values of local observables 
are independent of the boundary vectors Vr and Vr. 


As explained in m, to parametrize gauge-invariant 


MPS, i.e. states that obey G{n) |T(A)) = 0 for every 
n, it is convenient to give the virtual indices a multiple 
index structure a (q^aq); (3 {r,/3r), where q resp. r 

labels the eigenvalues of L{n — 1) res p. L {n). One can 
verify that the condition G{n) = 0 ( |2.5| ) then imposes 
the following form on the matrices: 


(2.9) 

where aq = 1... Dq{n), jSr = 1 ... D r{n + 1). The first 
Kronecker delta is Gauss’ law ( |2.5| ) on the virtual level 
while the second Kronecker delta connects the virtual in¬ 
dex r with the physical eigenvalue p of L{n). Because the 
indices q (resp. r) label the eigenvalues of L{n — 1) (resp. 
L{n)) and we only retain the eigenvalues of I/(n—1) in the 
interval Z[p^i^(n),p^oa,(n)] (resp. of L{n) in the inter¬ 
val Z[pmin{'n^l),Pmax{'n^l)]) wc havc that Dq{n) = 0 
for q > Pmaxi'^) and q < Pmin{'^)- The formal total bond 
dimension of this MPS is D(n) = ^ Dq(n), but 

notice that, as ( |2.9| ) takes a very specific form, the true 
variational freedom lies within the matrices as,p{n) G 

£^Dq(n)xDrin-\-l) 


Gauge invariance is, of course, also reflected in the 
Schmidt decomposition (|2.8|): for states of the form (2.9) 


the Schmidt values can be labeled with the same dou¬ 
ble index a (q^ aq). More specifically, the Schmidt 
decomposition (|2.8|) now reads (see Appendix A and G): 


Pmax (n+l) Dq(n-\-l) 

i^(^))= E E 

q=Pmin{n-\-l) Oiq = l 

( 2 . 10 ) 

As before, we observe that taking a finite bond dimen¬ 
sion Dq{n-\-l) corresponds to a truncation in the Schmidt 
spectrum, now of the charge sector q. The choice for the 
different bond dimensions Dq{n-\-l) in the different sim¬ 
ulations should then be such that the discarded Schmidt 
values for each charge sector are sufficiently small. For 
our simulations with zero background, a{n) = 0, in m 
we could take Dq = 0 for \q\ > 3, i.e. Pmin = ~3 and 
Pmax = +3. For the simulations with a nonzero back¬ 
ground field we find that for the same accuracy it suf¬ 
fices to consider eight g-sectors. But - n ot s urprisingly 
given the first term in the Hamiltonian (2.6) - we find 
the relevant eigenvalues sectors of L{n) to be centered 
around a dominant sector po that can be shifted away 
from po = 0 for some sites n. The largest Schmidt value 
in each g-sector decreases as we move farther away from 
g = Po- When |g— po| ^ 4 all the Schmidt values Xq^aqi'^) 
are sufficiently small and we can safely take Dq = 0, i.e. 
Pmax ^ Po + 4 and Pmin ^ Po “ 4. Wc refer to appendices 
\A\ and [^ for more details on the weight of the differ¬ 
ent sectors for the different simulations, see in particular 


figs. |12a| and |l2b| and figs. |2Qc| and [2Qd| for some explicit 
examples. 

From the Schmidt spectrum (2.10) one can extract dif¬ 
ferent measures for the entanglement. In the following we 
will always use the von Neumann entropy [83]. For the 
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half-chain cut at site n, to which we associate the position 
z = l/2)a in physical units, we then have: 

Xq,a^{n)log[Xq^a^{n)]. ( 2 . 11 ) 

q aq 


III. ASYMPTOTIC STRING TENSION 


We first study the large distance behavior of the po¬ 
tential as captured by the asymptotic string tension 
ctq = lim^^+oo This is the quantity that indi¬ 

cates whether the probe charges are asymptotically con¬ 
fined (ctq 7 ^ 0) or not (erg = 0). For the Schwinger model 
erg has been computed analytically in the strong-coupling 
expansion [771I81H86]. At the numerical front the most 
successful computation to date used finite-lattice scaling 
methods in a Hamiltonian formulation m- An advan¬ 
tage of our MPS simulations is that in contrast to m we 
can directly work in the thermodynamic limit {N oc), 
leaving only the x ^ oo interpolation to extract the con¬ 
tinuum results. The challenge of taking this continuum 
limit now lies in the diverging correlation length ^/a (in 
lattice units), as MPS simulations require larger bond 
dimensions for growing correlation length [53] . 


To find the asymptotic string tension we put a probe 
charge —gQ at — oo and a probe charge gQ at +oc. As we 
explain in the previous section, a probe charge pair trans¬ 
lates to a background electric field a{n) in the Hamilto¬ 
nian (2.6). In this case the background electric field will 
be uniform: a{n) = —Q^Mn. The Hamiltonian is then 
invariant under T^, a translation over two sites. In ac¬ 
cordance with this symmetry the appropriate MPS vari¬ 
ational ground-state ansatz takes the form 


|^(A(l),^(2))) 


Ed n A/^2n-l I 1^) 5 (^T) 

VnG^ / 


where 


l^n — i^^ntPn) ^ { I 5 1} ^ T ^\Pmaxij^ T 1)]? 

= |{«^n}n6z), vl G and 

Ai^{n) € £D{n)xD{n+i) ^a,kes the form (2.9) (n = 1,2). 
This corresponds to a general MPS (2.7) in the ther¬ 


modynamic limit {N +oc) where the tensors A^^{n) 
depend only on the parity of the site n: A/^ 2 n-i “ 1) = 

and A^^^{2n) = A^^^{2),\ln. As a conse- 
quence Dq{n),Pmin{i^) and Pmax{n) also depend on the 
parity of n. 

As we explain in appendix we were able to ac¬ 
curately approximate the ground state and its finite 
energy per site eg = fg/2A, with the total in¬ 
frared divergent energy, within the class of states 
Therefore, we perform imaginary time evolution {dr = 
idt) of the Schrodinger equation, idt |^(A(1), A(2))) = 



Q Q 


(a) (b) 



(c) (d) 

FIG. 1: (a): string tension erg. (b): electric field per 
site, (c): comparison with the strong-coupling result 
( |3.3| ) (dashed line) for m/g = 0.125 and m/g = 0.25 . 

(d): comparison with the weak-coupling result (3.4) 
(dashed line) for m/g = 1, 2,4. Inset: zooming in on the 
m/g = 4 curve. 


H |4^(A(1), A(2))), with the time-dependent variational 
principle (TDVP) [80l|82]. In appendix [A| we also explain 
how we chose the virtual dimensions {Dq{l), Dq(2)} and 
{Pmin/max{^)^Pmin/max{'^)} by investigating the Schmidt 
spectrum. In [57], we found the energy of the vacuum 
So = 2Neo for the zero-background field a{n) =0. In 
the same fashion, we now compute the string tension erg 
as the extra energy density, induced by the uniform back¬ 
ground electric field: erg = {Sq — So)/L where L is the 
length of our lattice. In units ^ = 1 we have L = 2N/ 
and therefore (7q{x) = \/x{eQ — cq). 

From the numerical point of view it is important to 
take the convergence of this UV-finite quantity (tq{x) 
as criterion for halting the imaginary TDVP time evo¬ 
lution. As we explain in more detail in appendix [^ 
we computed values for (7q{x) in this way, for x = 
100, 200, 300,400, 600,800 and perform a polynomial ex¬ 
trapolation in l/^/x similar to [79]. This indeed allows 
us to recover a finite value for hma,^oo thereby ex¬ 

plicitly verifying that the UV divergencies in the energy 
densities ^/xeq and ^/xeo cancel out. 

In fig. I^we plot our result for the continuum string 
tension erg computed for different values of the mass m/g 
as a function of the charge gQ of the external quark- 
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antiquark pair. Note that we only consider Q-values 
G [0,1[ as the string tension is periodic in Q: Q ^ Q —p 
upon L{n) L{n) + p for p G Z in the Hamiltonian 
( |2.6[ ). Note also that one can combine this transforma¬ 
tion for p = 1 with a CT-transformation {C = charge 
conjugation): 

L{in^ —y 1 — Ij{ti 1) y — 0{in -|- 1) 

(Jz{n) -(Jz{n + 1) cr^(n) ^ cr^(n + 1). (3.2) 


This transformation gives Q 1 — Q in the Hamil¬ 
tonian ( |2.6[ ) and therefore cjq = So for our 

calculations we can restrict ourselves to values Q G 
[0,1/2]. In practice we consider the explicit values: 
Q = 0.05, 0.10,0.15,..., 0.45,0.47, 0.48,0.5 and perform 
an interpolating fit. 

Our considered values ioi mjg interpolate between the 
strong- and weak-coupling regime. In the strong-coupling 
regime mjg the string tension is computed in mass 
perturbation theory from the bosonized field theory up 
to order 0((m/p)^) [85] 




m 


Z(1 — cos(27r(5)) 


nn?T!^E^'K 


4p2 


(1 — cos(47r(5)) 
(3.3) 

where Z = 0.15993, E+ = —8.9139. As one can observe 
in fig. Ic for m/p ^ 0 our results indeed converge to 
this analytic result that is plotted with a dashed line for 
m/g = 0.125 and for m/g = 0.25. 

In the weak-coupling regime g/m <C 1 we can easily 
compute the string tension in standar d pe rturbation the¬ 
ory from the continuum Lagrangian (2.1) (see appendix 
[b| . Up to order (p/m)^ we find the string tension for 

Q < 1/2: 


p^ 2 \ m? Gtt ) 


(3.4) 


with the value for Q > 1/2 following from the identifi¬ 
cation (jQ = (Ti-q for the compact formulation of QED 2 
that we are considering. In fig. [^one can again observe 
the convergence of our numerical results to this analytic 
result, now for p/m ^ 0. Not ice here that we subtract 
the leading order term of (3.4). 

Comparing the strong- and weak-coupling regime we 
observe an important difference: in the strong-coupling 
limit (jQ is differentiable at Q = 1/2 whereas in the weak- 
coupling limit this is not the case. Therefore there exists 
a critical mass (m/p)c with the property that cfq is differ¬ 
entiable at Q = 1/2 for {m/g) < {m/g)c and not differen¬ 
tiable at Q = 1/2 for (m/p) > {mjg)^^ This point {m/g)c 
corresponds to the first-order phase transition for the 
Hamiltonian Hq (2.6) at Q = 1/2 [7^. H q=ii2 is sym¬ 
metric under the Ci transformation ( |3.2[ ) and the point 
{m/g)c separates the unbroken phase m/p < {m/g)c 
from the spontaneously broken phase m/p > {m/g)c that 
was originally predicted by Coleman [44]. This relation¬ 
ship of the breaking of CT-symmetry with the nondiffer¬ 
entiability of erg can be made more concrete by noting 



FIG. 2: (a): m/g = 0.25, Q = 0.45. Fit of the form 
(-1/6) log(l/v^) + A + C/^/x to Sq{x) and So{x). 
Inset: linear fit to ASq{x). (b): ASq for different 
values of m/p. 


that 




\n=l,2 



(3.5) 


where (•. .)g denotes the expectation values with respect 
to the ground state of Hq. We now have the relation 
Eq = —Ei_q from the CT-transformation (3.2), which 
indeed makes it a good order parameter for the CT break¬ 
ing at Q = 1/2. 

We perform an independent computation of Fg, again 
for Q = 0.05, 0.10, 0.15,..., 0.45, 0.47, 0.48, 0.5, and now 
using values x = 100, 200, 300,400 for our continuum ex¬ 
trapolation (see appendix [a]). Our results are displayed in 
fig. lb At Q ^ 1/2 we find for m/p = 0.3, Eq/g = 0 up 
to a numerical error of 4 x 10“^ while for m/p = 0.35 
we find Eq/g = 0.314(2), consistent with the value 
{m/g)c ~ 0.33 that was obtained in 
sistent with the behavior of erg in fig. 


m 


la 


and also con- 


Fin ally we also compute the half-chain entropy S 
(2.11) for different values of Q and m/p, which in this 
translational-invariant case does not depend on the po¬ 
sition of the cut. As such the entropy is a UV divergent 
quantity, but one expects the divergence t o com e from the 
fermion kinetic term in the Hamiltonian (2.6) and there¬ 


fore be Q-independent. Specifically, the general results 
of Cardy and Calabrese [87] predict for two fermionic de¬ 
grees of freedom a UV divergence (with correlation length 
^ in physical units): 


1 


= — - log{l/^/x) -j- (finite terms as x 

6 


(3.6) 

+00) 


with a = Xjgy/x the lattice spacing. This is precisely 
what we find in our simulations. As an illustration, in fig. 
^we show a fit of the form (—1/6) log(l/v^)+A+C/v^ 


through our data of Sq{x) for Q = 0, Q = 0.45 and 
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mig = 0.25. A and B are obtained by a l inea r fit through 
Sq{x) + (1/6) log(l/v^), see subsection A3 of appendix 
[A| There, we also explicitly extract the coefficient —1/6 
of the logarithmic term by a logarithmic fit to the data. 
The errors are of the order 10“^ for m/ g > 0.5 and only 
of order 10“^ for m/g < 0.5, see table IV in appendix 


The universality of the logarithmic UV divergence then 
allows us to define a UV-finite renormalized entropy 
ASq = Sq — Sq^ with a finite continuum value that can 
be obtained by a polynomial extrapolation in 1/see 
inset fig. Contrary to the string tension and the elec¬ 
tric field, we found sometimes that the results at x = 100 
and the continuum results differ by a factor of order one 
or have different sign. We refer the reader to subsection 
|A3| in appendix and, in particular, to fig. for the 
details about the continuum extrapolation. In fig. [^we 
show this renormalized entropy ASq as a function of Q 
for different values of m/g. Most notably we observe an 
(almost) divergent behavior for m/g = 0.3 at Q 1/2 
close to the critical point Q = 1/2, {m/g\ ~ 0.33. From 
( |3.6[ ) we indeed expect a growing entropy for growing 
correlation length. By the same argument one can un¬ 
derstand the behavior at small Q-values: there the corre¬ 
lation length (inverse mass gap) increases with growing 
g/m [88], which is indeed paralleled by the behavior of 


IV. FROM SMALL TO LARGE DISTANCES 


Now, we consider the situation where the external quark 
and antiquark pair are separated over a finite length L. 
On a lattice with spacing a and interquark distance L = 
ka, the pair introduces a nonuniform background el ectri c 
field a{n) = —(30(0 < n < k) in the Hamiltonian (2.6). 
As ansatz for our MPS trial state |^(^)) for the ground 
state we now write: 


i$(B)) = ^q 

K 



Vr-1 \ 

n 

.n=rL / 




VR Ik) , 


(4.1) 


where tr < k tr and Ai^(n) =_^{n mod 2) cor- 
responds to the MPS approximation of the ground 
state of the zero-background Hamiltonian {a{n) = 0) and 
depends only on the parity of n. This is a MPS (2.7) 
in the thermodynamic limit {N Too) where we take 
A{n) = B{n) for ri < n < vr — 1 an d take the A{n) cor¬ 
responding to the ground state ( |3.1[ ) for a{n) = 0 to the 
left and to the right of the H(n)’s (n < vr and n > vr). 


The idea behind this ansatz is that the nonuni¬ 
form background electric field changes the vacuum and 
breaks translation invariance (all B{n) are different) 
but that asymptotically (|n| ^ 1) it does not af- 


VQ{L)/g (^(^)7°^'(^)) (Q = 1.5) 



(a) (b) 

FIG. 3: m/g = 0. (a): Potential for (3 = 1 and Q = 1.5 
compared with exact result in the c onti nuum ( |4.2[ ). 
Inset: convergence for x +oo to (4.2) for (3 = 1. (b): 
Distribution of fermion charge for Q = 1.5 for different 
separation lengths of the quark and antiquark for 
X = 100. The results are compared with the exact result 

(p). 


feet the vacuum. Again, gau ge in variance ( |2.5[ ) is im¬ 
posed if B{n) takes the form ( |2.9[ ) with general matri- 
ces 6s,p(n) e C°-i(n)y<D^(n+l) g ’Z[pmin{n),Pmax{n)]; 
p,r e Z\pmin{n + + I)])- Note that we allow 

different bond dime nsions on different sites. Also as^p{n) 
is of the form ( |2.9| ) as we impose this to determine the 
ground state of the zero-background electric field Hamil¬ 
tonian. 


Because (4.1) is linear in each of the B^ we can use 
the DMRG-method [89] to obtain the best approxima¬ 
tion for the ground state within the manifold of gauge- 
invariant states, by optimizing on the UV- and IR- 
finite quantity Vq{L). By looking at the Schmidt spec¬ 
trum we are able to fix the values of the virtual dimen¬ 
sion Dq(n) and the minimum and maximum eigenval¬ 
ues Pmin{n + 1) and Pmax{'^ + 1) of I/(n) wc retain in 
our numerical scheme to obtain an accurate approxima¬ 
tion of the ground state. The choices for tr and vr, 
which vary between —k/2 — 250 < vr < —k/2 — 150 and 
k/2 -\- 150 < vr < k/2 250 are checked a posteriori 

by verifying the convergence of local observables at large 
distances to their zero-background value. We refer the 
reader to Appendix [C] for the details. 


A. The case m/g=0: screening a la Higgs 


We first discuss our results for the m/g = 0 case. This is a 
special case, as the asymptotic string tension ag vanishes 
for all values (integer or fractional) of the charge. Phys¬ 
ically, this is interpreted as a manifestation of a Higgs 
mechanism suppressing the long-range Coulomb 

force and replacing it with a short-range Yukawa force 
thereby effectively screening all charges. Another reason 
that makes the m/g = 0 case special is that it can be 
solved analytically [84], which allows for benchmarking 




















0.03 


ASq{z} {L-g = 0.85) 


0.04 


A5q(z) {L-g = 5.25) 




(a) (b) 

FIG. 4: m/g = 0, Q = 1. Spatial profile of ASq for 
different values of L and scaling to the continuum limit 
(x +oc). (a) Lg = 0.85. (b): Lg = 15.65. 


of numerical results. Previous numerical calculations for 
this case were performed with Monte Carlo simulations 
on the bosonized version of the theory m- 

In fig. we plot our results for the potential for 
m/g = 0 for Q = 1 and Q = 1.5. This can be com¬ 
pared with the exact continuum result [84]: 

(L) = (^1 _ e-^s/, (4.2) 


which is indeed of the Yukawa type. We find very good 
agreement already for x = 100, for both Q = 1 and 
Q = 1.5. For Q = 1 we also perform a computation for 
X = 400, in the inset of fig. [^one can observe the rate of 
convergence towards the continuum x ^ oc in this case. 

The charge density of the light quarks 

is, of course, also an interesting quantity to compute, as 
it explicitly shows the screening of the external probe 
charges. The analytical result for the probe charge pair 
put at ±L/2 reads [84] : 


{^p{z)'y°^p{z)) = 



^Q-9\z+L/2\/y/^ _ ^-g\z-L/2\/y^^ ^ 


(4.3) 

This indeed corresponds to a charge distribution with two 
‘clouds’ of oppositely charged light (in this case massless) 
quarks, around the external quark and antiquark, that for 
large distance L have exactly the same total charge AQ 
as the external pair. On the lattice the charge density at 
2 ) = (2n—l/2)a is computed as ^/x{az{2n—l)^az{2n))/2. 
In fig. we plot this density for Q = 1.5 where the 
charges are separated at distances Lg = 5.1 and Lg = 
17.3. Here, too, our results for x = 100 are already very 
close to the continuum result. 


In fig. we show the spatial profile of the renor¬ 
malized half-chain von Neumann entropy ASq{z) = 
Sq{z) — So{z) for different values of Lg. We compute 
this quantity for z = {n-\- l/2)a with n even and perform 
an interpolating fit. When the heavy quarks are close 
to each other, ASq{z) shows a peak in the middle be¬ 
tween the charges and falls off very fast with \zg\^ see fig. 


4a For larger values of a cloud of light quarks forms 


around each of the heavy charges which clearly leaves its 
imprints on the spatial profile of the von Neumann en¬ 
tropy, see fig. |4b[ ASq{z) is nonzero around each of the 
heavy charges and is zero around zg ^ 0. 

The observed spatial profiles of the von Neumann en¬ 
tropy, however, are lattice artifacts and vanish in the 
continuum limit {x +oo). Indeed, from the bosonized 
Hamiltonian for m/g = 0, it can be observed that 
any position-dependent electric background field can be 
transformed away na E] up to a position-dependent 
constant. Therefore the von Neumann entropy of the 
ground state with a{n) ^ 0 and a{n) = 0 are the same, 
hence ASq{z) = 0 for any value of Lg. By investigating 
the scaling towards x +oo in figs. [4a|and[4l^we indeed 
observe that ASq{z) tends towards a very small value for 
X Too. Note that we here need to perform an interpo¬ 
lation because we can only take to be an integer mul¬ 
tiple of l/\/x. Specifically, we perform simulations for 
X = 400 and Lg = 0.85, 5.25,15.65. For x = 100, 200, 300 
we first do simulations for Lig < 0.85, 5.25,15.65 and 
1/2^ > 0.85, 5.25,15.65. Afterwards we do a simple linear 
interpolation between Lig and L 2 g to obtain the curve 
for Lg = 0.85, 5.25,15.65. By performing a linear ex¬ 
trapolation in 1/v^ through x = 100,200,300,400 and 
X = 200, 300,400 we find that in fig. |^the continuum ex¬ 
trapolation of the maxima for zg ^ 0 yields the estimate 
ASo(O) ~ 5(7) X10“^, while extrapolating the maxima in 
fig. 4b around zg = ±10 gives A5 'q(± 10) « 2(5) x 10“^, 
consistent with ASq{z) = 0. The interpolation to ob¬ 
tain the curves for Lg = 0.85, 5.25,15.65 for all values of 
X leads to relative large errors such that a continuum ex¬ 
trapolation of ASq{z) for all values of zg is less reliable. 
Therefore, for Lg < 0, we perform a linear extrapolation 
in 1/v^ through our values for x = 100, 200, 300,400 and 
obtained the results for > 0 by a reflection of the re¬ 
sults for Lg < 0 (green line). This can be compared with 
the analytical result ASq{z) = 0 (orange). 


B. The case Q = 1: string breaking 


For m/g ^ 0, the asymptotic string tension cfq vanishes 
only for integer charges Q. This is taken to be an indi¬ 
cation for a screening a la QCD where the potential 
exhibits a string tension (Vq{L) oc L) at short distances, 
but flattens out completely at large distances, at least for 
integer charges Q. At these large distances it becomes 
energetically favorable to materialize light (yet massive) 
(anti)fermions out of the vacuum that bind to the exter¬ 
nal quark and antiquark, resulting in two charge neutral 
mesons. 

Historically, for QCD, lattice Monte-Carlo simulations 
succeeded first to calculate numerically the short distance 
confining behavior of the potential - both in the quenched 
and unquenched approximation - via the expectation 
value of the Wilson mm- The detection of string break¬ 
ing has posed a larger challenge. A main problem with 
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(a) (b) 



(c) 


(d) 



(e) (f) 

FIG. 5: Q = l,x = 100 . (a): Quark-antiquark potential 
for different values of m/^. (b): Comparison of potential 
with nonrelativistic limit result (dashed line) for m/g = 
0.25, 0.75, 1, 2. (c): The total charge of the light 
fermions on the negative axis Q-{L) for different values 
oimjg. (d): Electric field for m/^ = 0.75. (e): Charge 
distribution for mjg = 0.75. For Lg = 17.3 we compare 
with the charge distribution of the nonrelativistic meson 
state (full red line), (f): Comparison of the charge 
density of the left cloud (full li ne) w ith that of the 
nonrelativistic meson state, eq (4.5), (dashed line) for 
Lg = 17.3, now for mjg = 0.125, 0.5,1, 2. 


the use of the standard Wilson loop is the poor overlap 
with the broken-string two-meson state. This problem 
was finally overcome by including light quark propaga¬ 
tors in the Wilson loop and analyzing its mixing with 
the standard Wilson loop [92] . 

For the Schwinger model the string-breaking phe¬ 


nomenon has been confirmed in mass perturbation theory 
j86] and in a semiclassical approximation of the bosonized 
version of the theory [93l|94]. At the numerical level, for 
(5 = 1, lattice Monte Carlo simulations have detected 
both the confining and string-breaking behavior of the 
potential [90l [95]. In [90] the problem with the Wilson 
loop was avoided by computing instead the expectation 
value of the bosonized Hamiltonian, while [95] turned to 
very high statistics thereby explicitly showing the poor 
overlap of the Wilson loop with the broken-string ground 
state. 

For the local quantities (charge density =7/5(^) 7 ^ 7 /^( 2 :), 
electric field=F^( 2 :)) and the potential, we restrict our¬ 
selves from now on to lattice spacing x = 100(= Xjg^o?)'^ 
from the previous subsection we can expect these results 
already to be quite close to the continuum. In fig. |5^we 
display our results for the potential, and this for different 
values of mjg. We compute explicitly the ground-state 
energy at = 0.1,0.3,... 15.3 and perform an interpo¬ 
lating fit. We clearly find a transition from the confining 
behavior, associated with the string state, towards the 
constant behavior associated with the broken-string two- 
meson state. This transition happens more suddenly for 
larger values of m/^, which is in qualitative agreement 
with the semiclassical results from the bosonized theory 
|93l [94]. This is also what one would expect from the 
nonrelativistic weak-coupling regime, where the transi¬ 
tion can be understood as a level crossing between the 
zero-particle string state and the two-particle broken- 
string meson state (see appendix 0. The dashed lines 
in fig. corresponding to this nonrelativistic result for 

^string = Lg^ 1‘i and £2meson = 2 m + 1. 0188^173 were 
plotted for comparison. We can indeed observe the con¬ 
vergence towards this result for increasing values oimjg. 

We further illustrate this behavior in fig. where 
we plot the total charge Q- of the light fermions on the 
negative axis: 


Q- = 9 



dz . 


(4.4) 


One can observe indeed that the interpolation between 
(5- = 0 (string state) for small L and Q- = 1 (meson 
state) for large L becomes more and more discontinuous 
for growing mjg in accordance with the nonrelativistic 
level-crossing picture. 


In figs. [5d|and fig. [^we investigate the interpolation 
from the string state to the string-broken state in more 
detail for mjg = 0.75 by plotting the charge density and 
electric field. For Ljg = 0.5 there is only a very small 
charge cloud around the external quark and antiquark, 
notice also the very short electric field string displayed 
at the bottom of fig. 5d At Ljg = 5.1 the clouds start 
to build up, lowering the electric field value at the cen¬ 
ter. Ai Ljg = 10.1 the string is completely broken, the 
electric field at the center has vanished, and we have two 
clouds of total charge ±1 around the external quark and 
antiquark. At Ljg = 17.3 the two isolated mesons are 
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(a) (b) 

FIG. 6: m/g = 2^ Q = 1. ASq{z) for different values of 
L and scaling to the continuum limit, (a) Lg = 5.25. 

(b) Lg = 10.95. 

simply separated over a larger distance, with a quasi¬ 
identical charge distribution around the external quarks 
as for Ljg = 10.1. 

The full red line in fig. is the charge distribution 
for the nonrelativistic meson state for Lg = 
17.3, with 

(l){z) = AfAi (^{g^my^^\z ± L/2\ — 1.0188^ (4.5) 

the ground state of the one-particle problem in a linear 
potential (appendix 0 , where Ai is the Airy function 
m and JV the normalization factor. As one can observe, 
the charge distribution from this nonrelativistic picture 
matches very well our exact (numerical) result. In fig. 
we compare the charge cloud at the negative 2 ;—axis with 
the nonrelativistic result for other values of m/g. One 
can again observe the convergence to the nonrelativistic 
result for growing m/g, notice that already for m/g = 0.5 
the match is quite good. 

For the renormalized von Neumann entropy ASq{z), 
we also find a characteristic picture, both for the string 
state and the string-broken state, see fig. [^for the case 
m/g = 2. For the string state, Lg < 9.5, the entropy 
shows a constant surplus in between the probe charges, 
similar to the electric field. But notice that this effect be¬ 
comes very small in the continuum limit (green line), we 
find an extrapolated value: ASg{z) ~ 2.0(5) x 10“^ for 
zg G [—2.5, 2.5]. For the string-broken case Lg > 10, see 
fig. we find that the entropy now shows two clouds 
around the heavy quark and the heavy antiquark, sim¬ 
ilar to the charge density. But notice that in contrast 
to the string state the entropy now survives the contin¬ 
uum limit, with the x = 100 value already close to the 
continuum extrapolation. 


C. General Q: partial string breaking 


We now finally turn our attention to the general case 
Q 1. In this case we should have the interesting phe- 



(a) (b) 

FIG. 7: Gartoon picture of string breaking in the 
nonrelativistic limit m/g ^ 00 . The electric field gets 
successively screened by light quarks (antiquarks) that 
bind to the external charges with charge -bQ (-Q). This 
leads to the formation of two mesons: one meson 
existing of the heavy quark and light antiquarks with 
charge +1 and one meson existing of the heavy 
antiquark and the light quarks with charge -1. (a) 

Q = 4.5: the remaining electric field in between with 
net charge —0.5 confines the two meson configurations 
asymptotically, (b) (5 = 5: the electric field is entirely 
screened and the meson configurations are deconfined. 


nomenon of partial string breaking. Indeed, in the non¬ 
relativistic limit m/g 00 of string breaking due to 
meson formation, probe charges Q can only be screened 
by an integer number: Q ^ Q = Q — n, where n is 
the number of light (anti-)quarks that bind to the exter¬ 
nal charges. For nonzero Q, i.e. when Q is noninteger, 
this still leaves a string between the two separated me¬ 
son configurations. A visualization of this process in the 
nonrelativistic limit is shown in fig. 

Our simulations allow us to verify to what extent this 
picture is realized for finite m/g. In fig. |^we plot our re¬ 
sults for different values of Q, both fractional and integer. 
We do indeed recover partial string breaking, largely fol¬ 
lowing the nonrelativistic picture. To our knowledge this 
is the first successful simulation of partial string breaking 
in the Schwinger model, a previous Monte Garlo simula¬ 
tion [ 95 ] failed to detect the phenomenon. 


In fig. [^we plot, as in the previous section, the evo¬ 
lution of the total dynamical charge Q- at the negative 
z-axis, for m/g = 1. For all values of Q this charge Q- 
indeed makes quasidiscrete jumps of AQ_ ^ +1 which 
should correspond to (partial) string breakings. As we 
see in fig. these jumps indeed correlate with jumps in 
the string tensions in the different regions of the poten¬ 
tials. For m/g = 0.5 we still find jumps of Q- but they 
are smoothened out, as can be seen in figs. and |8d| 
For m/g = 0.25 the jumps are even more smoothened 
out as can be seen in figs. and This smoothened 
behavior, similar to what we obtain in the Q = 1 case, is 
expected as we go further from the nonrelativistic large 
m/g regime. But still note the contrast with the be havior 
in the massless limit m/g = 0 of subsection IV A, where 
the charge Q- grows continuously to the external value 
Q, assuring a complete screening. 

For L going from 0 to oc, different partial string break- 
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VQ(L)lg (m/<,= l) 




(a) (b) 



(c) (d) 



(e) (f) 

FIG. 8: X = 100. (a): Quark-antiquark potential for 
m/g = 1 for different values of Q. (b): Q-(L) for 
m/g = 1 for different values of Q. (c) and (d): the same 
quantities for m/g = 0.5. (e) and (f): the same 
quantities for m/g = 0.25 


ings should lead to the asymptotic behavior of the poten¬ 
tial that we examine in section Hill In table [T] we show the 
difference of the slope of the potential around Lg = 15.3 
with the asymptotic string tension at x = 100 that we 
calculate in the previous section. The former is estimated 
as the mean of the backward differences 

1 AFq _ VQ(15.3g) - VQ(Lg) / 1 cIVq) 

AL Lg r 9 ^ dL ) ^ 

for Lg = 13.7,14.1,14.5,14.9. The error is computed 
as the standard deviation of these backward differences. 
One observes that for m/g = 1 the string tension has 



m/g 

Q 

0.25 

0.5 

1 

0.75 

2.5(6) X 10“® 

5(2) X 10“’' 

-7(2) X 10“^° 

1 

8(2) X 10“’' 

5(2) X 10“® 

-4(4) X 10““ 

1.75 

2.6(7) X 

1.7(7) X 10“® 

8(8) X 10“® 

2.5 

3.0(1) X 10“® 

2.5(5) X 10“® 

-1(1) X 10“*^ 

3.25 

2.2(1) X 10"® 

2.1(8) X 10“’’ 

-1(1) X 10“*^ 

4.5 

4.0(2) X 10“® 

1.0(2) X 10“® 

-1(1) X 10“® 

5 

2.1(6) X 10"^ 

1.0(5) X 10“® 

-2(2) X 10“® 


TABLE I: X = 100. Values for the difference 
{AVq/AL — (yQ)/g^ where AVq/AL is the mean of the 
backward differences at Lg = 15.3 with AL g = 
0.4,0.8,1.2,1.6. 



m/g 

Q 

0.25 

0.5 

1 

0.75 

0.9225 

0.9675 

0.9891 

1 

0.9995 

1.0000 

1.0000 

1.75 

1.9157 

1.9665 

1.9891 

2.5 

2.2384 

2.0778 

2.0223 

3.25 

3.0748 

3.0332 

3.0111 

4.5 

4.2150 

4.0770 

4.0230 

5 

4.9922 

4.9990 

5.0000 


TABLE II: x = 100. Values for Q- at Lg = 15.3 for 
m/g = 0.25, m/g = 0.5 and m/g = 1 . 


already converged to the asymptotic result, almost up 
to the numerical precision, while for m/g = 0.5 we 
are already very close to the asymptotic result and for 
m/g = 0.25 there is a slightly larger (but still very small) 
difference. 

Eor integer values of Q, the asymptotic string tension 
vanishes, so asymptotically we expect Q- ^ Q, corre¬ 
sponding to a complete screening. Eor the values Q = 1 
and Q = 5 that we consider, this is already almost sat¬ 
isfied at Lg = 15.3, as can be seen in table [Ilj In the 
nonrelativistic limit for general Q, the total dynamical 
charge Q- that is produced asymptotically, will be the 
integer number that minimizes \Q — Q-\. Eor finite m/g 
we expect corrections to the nonrelativistic limit, but as 
one can see in the table these corrections are still very 
small for m/g = 1 and m/g = 0.5. Notice also that for 
the half-integer values (3 = 2.5 and 4.5, for which we have 
spontaneous symmetry breaking in the asymptotic limit 
(see section [lll| ) , we find Q- approaching the smallest of 
the two possible nonrelativistic values Q- ^ Q — 1/2. 

In fig. [^we show the spatial charge distribution and 
electric field for different distances of the probe quarks. 
Eor Q = 1.75 we have two partial string breakings. 
The first one, around Lg « 1.7 (see fig. brings the 
electric field string at the center from E/g ^ —1.7 to 
E/g ^ —0.7. After the second partial string break¬ 
ing, around Lg ^ 9, the probe charge is ‘overscreened’, 
(3- ~ 2, leading to a final electric field string with op¬ 
posite sign E/g ^ -1-0.2. Notice that in contrast to the 
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(a) (b) 



(c) (d) 



(e) (f) 

FIG. 9: mjg = 0.5, x = 100. Left: Charge distribution 
for Q = 1.75 (a), Q = 4.5 (c) and Q = 5 (e) for different 
values of the separation length L. Right: Electric field 
for Q = 1.75 (b), Q = 4.5 (d) and Q = 5 (f) for different 
values of the separation length L. 


Q = 1 case, the charge clouds at large separation of the 
probe quarks are not symmetric around the position of 
the probes. This is expected, as the remaining confin¬ 
ing force between the two (charged) ‘mesons’ distorts the 
charge distribution. For Q = 4.5 we have a similar pic¬ 
ture, but now, after the final partial string breaking, the 
probe charge is ‘underscreened’, Q- ~ 4, resulting in a 
final negative electric field string E/g ^ —0.4. While for 
(5 = 5 the final string breaking is complete: the probe 
charge is screened entirely (5- ~ 5, leading to a complete 
neutralization of the electric field string E/g ^ 0 at the 
center. In this case for large enough Lg we expect the 


ASq{z) {L-g = 0.55) ASq{z) {L ■ g = 2.55) 




(c) (d) 

FIG. 10: m/g = 0.5, Q = 4.5. ASq{z) for different 
values of L. We also show the scaling to x ^ -hoc. (a) 
Lg = 0.55. (b) Lg = 2.55. (c) Lg = 7.35 (d) Lg = 13.25. 


charge distributions to become fully symmetric around 
the probe charge positions. 

In fig. [^we show the effect of different partial string 
breakings on the entropy profile ASq{z), for m/g = 0.5 
and Q = 4.5. For the smallest interquark distance 
Lg = 0.55, the entropy peaks at the center. At Lg > 2.55 
(after two string breakings, see fig. [^, we observe a pro¬ 
file with two peaks around the positions of the probe 
charges. At Lg > 7.35 and Lg > 13.15, after four string 
breakings, the profile now shows four peaks around the 
probe quark positions. In addition, we find an entropy 
surplus in the center, which now seems to be stable under 
the continuum extrapolation. 

In fig. [^we show that this characteristic imprint on 
the entropy is generic. We plot ASq{z) for Lg = 15.25 
and different values of Q. For Lg = 15.25 all the par¬ 
tial string breakings have occurred and the final meson 
configurations around the external charge positions are 
formed. By counting the peaks one can again deduce the 
number of light elementary quarks (corresponding to the 
number of partial string breakings) in the meson states. 
For instance, for Q = 4.5, fig. |lla| we observe that 
there are four partial string breakings and for (5 = 5, fig. 
[iibl we observe that there are five partial string break¬ 
ings. The spatial profiles do in fact differ by only one 
additional peak in each of the clouds for (5 = 5 around 
zg = ±5. Notice also the difference in the spatial profile 
for Q = 1.75 and Q = 2.5, see figs. |llc| and |lld[ In both 
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(c) (d) 

FIG. 11: mjg = 0.5, Lg = 15.25. ASq{z) for different 
values of Q and scaling to x ^ +oc. (a) Q = 4.5. (b) 
Q = 5. (c) Q = 1.75. (d) Q = 2.5. 


cases two partial string breakings lead to the asymptotic 
meson state, but in the former case the final electric field 
is ‘overscreened’ while in the latter case the final electric 
field is ‘underscreened’. Finally, notice that we can trust 
these results to be close to their continuum value, as the 
variation for the different x-values is very small. 


V. CONCLUSIONS 


In this paper we employed the MPS formalism for a de¬ 
tailed numerical study of the confining mechanism in the 
static limit for the massive Schwinger model. Our Hamil¬ 
tonian setup gives us direct access to the modified vac¬ 
uum state in presence of two probe charges. This allowed 
us, not only to compute the interquark potential, but 
also the spatial profile of the electric field between the 
probe charges and the charge concentration of the light 
fermions. Even for relatively small m/g the picture that 
emerged can be understood as a smoothened version of 
the nonrelativistic limit, with a level crossing between 
the electric string state that is the ground state at short 
distances and the broken-string two meson state that is 
the ground state at large distances. Here the two isolated 
mesons each consist of a light (anti-)quark cloud around 
the heavy probe charge, that is well described by the 
solution to the Schrodinger equation of the appropriate 
one-particle problem. 


In the case of fractional probe charges, we clearly 
observed the expected partial string breaking. Again 
in accordance with the nonrelativistic picture we found 
the screening of the probe charges to happen in jumps 
AQ ^ 1 of the light fermion charge; with these jumps 
becoming more and more discrete for growing m/g. 

Our tensor network simulations also give us direct ac¬ 
cess to the full Schmidt spectrum for the different biparti¬ 
tions on the state. The numerical simulations show that 
the UV divergence in the corresponding von Neumann 
entropy is universal, allowing us to define a UV-finite 
renormalized entropy by subtracting the vacuum value. 
We have examined the imprint of both the string forma¬ 
tion and string breaking on the profile of this renormal¬ 
ized entropy. Most notably we found that string breaking 
leaves a very distinct imprint on this entropy profile. 

We have checked our results not only against the 
predi ctions from the one-particle Schrodinger equation 
(D.3), but also against the weak-conpling results from 


the original Lagrangian (2.1) and against the strong- 
coupling results from the bosonized field theory [44]. In 
the appropriate regimes we found nearly perfect agree¬ 
ment with these continuum analytic results. This not 
only demonstrates the potential of MPS simulations close 
to the continuum critical point of a lattice theory. But it 
also serves as a nice, if not unexpected, cross-check of the 
consistency of all different descriptions of the Schwinger 
model. 


We have restricted ourselves to the study of the static 
limit of the confinement mechanism. An obvious future 
extension of our work is to consider the dynamical prob¬ 
lem, simulating the real-time hadronization that takes 
place in a realistic scattering process. MPS real-time 
simulations of this type of problem were considered re¬ 
cently for U{1) and SU{2) quantum link lattice models 
in [58] and m- One could also approach this problem by 
first calculating the scattering eigenstates [97l [98] . See 
also (99] llQQj for an approach in the semiclassical limit. 

Of course it will also be very interesting to bring 
this type of analysis to higher dimensions. Specifically, 
present techniques should already allow one to simu¬ 
late e.g., the static confinement problem for some simple 
2+1 dimensional lattice models (see m for homogeneous 
ground-state simulations of a Z 2 model and m for simu¬ 
lations of a 1/(1) lattice model with dynamical fermions). 
However, the current algorithms for PEPS simulations 
scale unfavorably with the bond dimension and we 
therefore expect that the successful simulation of specific 
microscopic gauge field Hamiltonians in the continuum 
limit will require new techniques. Still, given the contin¬ 
uous progress of PEPS methods |1Q21[1Q5] we are hopeful 
on that front. In any case, in light of the potential for 
real-time and finite fermion density simulations, and of 
the new insight that might come from understanding the 
entanglement structure, it should certainly be worthwhile 
to further explore this direction, hopefully succeeding one 
day in the full simulation of the microscopic Hamiltonian 
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of (3+l)-dimensional QCD. 
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Appendix A: Computation of the asymptotic string tension and electric field 


In this appendix we d iscuss the details of the computation of the asymptotic string tension erg and the electric field 
E. In subsection A 1 we discuss how to ob tain the electric field E{x) and the energy density e{x) for a fixed lattice 
spacing ga = Ij^Jx and in subsection A 2 we discuss how we to extrapolate these quantities to the continuum limit 
{x +oc). 


1. MPS ansatz 


To find the electric field E{x) and the energy density t{x) for a fixed lattice spacing ga — Ijy/x we need to find the 
ground state of the Schwinger Hamiltonian 


= A ( + (-1)") + xY,{a+{n)e^^^-^a-{n + 1) + ft.c.)) (A.l) 

^ VnG^ nG^ nG^ ^ 


in an electric background field Q G [0,1[. Taking into account the translation symmetry over two sites, we propose 
the following MPS ansatz in the thermodynamic limit, see eq. (3.1), 


|^(A(1),A(2))) = ^v\ ( ) Vi? |/«), 


(A.2) 


VnGZ 


with 


^ { I 5 1 }iPn ^ '^\Pminij^ T ^\Pmax{E ^)]5 1^) — |{^n}nGz) • 

Gauss’ law, G{n) = L(n) — L(n — 1 ) — {(Jz{n) + (— 1 )’^)/ 2 , imposes the following form for A(l) and A( 2 ), see eq. ( |2.9[ ), 

~ = 1 , 2 , (^-4) 

q. ^ ^[Prnin{n),Prnax{ri)]]p,r G Z[pmin{nEl) ,Pmax{nEl)],s = ±1, aq = l...Dq{n), pr = 1 ... (^ + 1 ) • The Optimal 

approximation for the ground state within the class of MPS with fixed bond dimension is obtained by performing 
imaginary time evolution of the Schrodinger equation using the time-dependent variational principle (TDVP) [82]. 
When applying the TDVP, the sites 2n — 1 and 2n are blocked into one effective site n\ 

(^S1,P1 (1): ^S2,P2 (2)) ^ ^S1,P1,S2,P2 ~ (1)^S2,P2 (2) (^•^) 

where 


[^Sl,pl,S2,P2]{q,CXq),{r,^r) ~ [^Sl,Pl,S2]aq,/3r^P2,Pl+(s2 + l)/2^P2,r^Pl,g+(si-l)/2 (^•^) 

with as^^p^^s2 ^ ]J^ = Dq{l). Note that = 0 for < pmin and q > Pmax where Pmin = Pmin(l) and 

Pmax = Pmax{^)‘ GaUSs’ laW implies that Pmin{‘^) ^ {Pmin - l,Pmin} and Pmax{‘^) G {p max ^iPmax} • 


To check whether the obtained MPS is a good approximation we should look at the Schmidt values Xq^aqi'^) 
associated with the bipartition {Mi(2n) = (Z[—oo, 2 n],M 2 ( 2 n) = Z[2n + 1, +oc]} of the lattice. Transla tion s ymmetry 
implies that the Schmidt values are independent of the effective site n: = Ag^Q,g,Vn. Similar to ( 2 . 10 ) we have 

then 


|^(A))= X] E (A.7) 

Q—Pmin <Aq = l 

where (resp. are orthonormal unit vectors in the tensor product of the local Hilbert spaces in the 

region Mi(2n) (resp. M 2 ( 2 n)). The Schmidt values Xq^^q^ which are non-negative and sum to one, can be obtained as 
follows: assume is brought in a canonical form such that the matrices r and I corresponding to the right and 

left eigenvectors of the largest eigenvalue of the transfer matrix [80] . 

( Pmax^k^ \ 

^ ^ ~ ^ ^ ^ ^ 5 ~ {^k:Pk) I ('A--8) 

Kk Sk = -iPpk=Pminik) J 
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are positive definite and diagonal. He re w e assume that the largest eigenvalue of the transfer matrix is normalized to 
one. Then, because A takes the form (A.6), r and I will also be degenerate in the eigenvalues of L{2n)\ [r](^q^a^)-(^r,^r) ~ 

The Schmidt values Xq^cxq are now obtained by multiplying r and /: 


^q,aq ^q,r^( 


'^q,aq^q,r^aq,/3r 5 [^] {q,aq)]{r,/3r) 
\,cxq = rq^aqlq,aq whcrC q G 
freedom of the matrices 


^Pmax] labels the eigenvalues of L(2n) 


and aq = 1... Dq labels the variational 


As can be observed from eq. (A.7), truncating to a finite bond dimension thus corresponds to an effective truncation 
in the Schmidt decomposition of the ground state. Ideally one would want a distribution of Dq-valnes such that the 
smallest retained Schmidt value is more or less equal for each eigenvalue sector of L(2n). Then if we want a reliable 
MPS approximation for the ground state, these smallest retained Schmidt values should be sufficiently small, which 
corresponds to taking Dq sufficiently large. In practice we do several simulations and adapted Dq until the smallest 


Schmidt value in each eigenvalue sector of L{2n) was of order 10 
In figs 


-17 


i.e. mmo 


A, 


10 


-17 


12a and 12b we plot the distribution of the Schmidt values among the eigenvalue sectors of L(2n) for the 


final MPS ground-state approximations for m/g = 0.75,x = 400 and Q = 0.2,0.45. As in [57], we observe that the 
sectors corresponding to q = 0, — I, I are the most dominant ones which justifies our ch oice of taking Dq = 0 for 
\q\ > 3. This can be understood physically from the term proportional to [L{n) — Q]‘^ in (A.I) which punishes l arge 


I2c 


eigenvalues of L{n). We also display the bond dimensions for each sector and for each simulated value of x in figs. 
and|I2d One can observe that as x increases we need larger Dq for the same accuracy. This is explained by the fact 
that the correlation length diverges as we approach the continuum limit {x Too) and it is well known that critical 
theories require larger bond dimensions for a good MPS approximation. For the same reason we also need larger Dq 
when we are getting closer to the phase transition at m/g = {m/g)c ~ 0.33 and Q = 1/2 [44l[79] . 


2. Continuum extrapolation of the string tension and the electric field 


In the second part of this appendix we discuss how we obtain an estimate for the continuum value of erg and Eq. 
Note that the string tension at x = Xjg^o? is obtained from the energy density by: 


oq{x) = Vx{eQ{x) - eo{x)) 


(A.9) 


where eg(x) is the ground-state energy per site of the Schwinger Hamiltonian (A.I). 
reduces to the XT-model we have that 


As for X ^ oc, H/{2g^/x) 


li.„ lim 


eo(x) 


-1 

TT 


a:^+oo 2g^/x cc^+oo 2g^/x 
and it is argued in m that eQ{x)/^/x should behave polynomially as a function of l/^/x for large x, we have 


(A.IO) 


V^"-^ = -- + CqV^ + Aq + o(2-] (x»i), 
9 TT ^ ^ V\/^/ 


(A.lla) 


9 TT \y/x; 


(A.llb) 


This means that the energy densities y/xeQ{x) and ^eo{x) are UV divergent. But as we see, the string tension 
which is the difference of these quantities is UV-finite a nd, thu s, we should also have Cq = Cq. However, from the 
numerical point of view it is clear that small errors in (|A.I0) or/and in Cq and Cq would lead to large errors in 
the extrapola ted c ontinuum value liuix^oo To avoid this problem we first calculate eo and subtract it from the 
Hamiltonian (A.I): Hq ^ Hq — The string tension is the n com pute as (tq{x) = gy/xeQ{x) where eg(x) is 

the ground state of the renormalized Hamiltonian. As follows from (A.II), for large x, cfq{x) should scale as 




\fx 


2,5/2 


(A.12) 


In our simulations we computed (Tq{x) for x = 100,200,300,400,600,800. Our estimate (Tg* is obtained by fitting 
the (jQ {x) corresponding to the five largest x to 


h{x) = Aq + ^ ^ 

yjx X 


Pq. 

2.3/2 


(A.13) 

























20 


logio \,a, {Q = 0-2, X = 400) 


logio Ka, iQ = 0-45, X = 400) 



(a) 

Dq [Q = 0.2) 




(b) 

Dq [Q = 0.45) 



(c) (d) 

FIG. 12: m/g = 0.75. (a): Q = 0.2, x = 400. Distribution of the 10-base logarithm of the Schmidt values 
among the eigenvalue sectors q of L(2n). (b): Same as (a) but now for Q = 0.45. (c): Q = 0.2. Distribution of the 
bond dimension among the eigenvalue sectors of L(2n) for different values of x. (d): Same as (c) but now for 
Q = 0.45. 














21 






FIG. 13: Q = 0.3 : Continuum extrapolation of the string tension erg for different values of m/g. 


and taking ag = q^Aq. 

In fig. 13 we plot our results for the string tension as a function of 1/^/x^oTQ = 0.3 and m/g = 0.125,0.3,0.5,1. 
The numerical results are represented by circles and our polynomial fit (A. 13) through the largest five x-values is 
shown by a full line. The star represents our continuum estimate. It is clear that the string tension indeed behaves 
polynomially as a function of l/\/x. For larger values of m/g one can also deduce that we are already very close to 
the continuum limit at l/^/x = 0.1. Indeed, for m/^ = 1, the difference of our estimate with (Jq{x) at x = 100 is only 
of order 10“^. 

The continuum extrapolation depends on the chosen interval and the chosen fit. Therefore we also compute the 
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(a) (b) 

FIG. 14: (a): log;LO (^^^ctq) ^ function of Q. (b): log;L 0 ^ function of Q. 


continuum estimates by fitting all the data to fi{x) (see ( |A.13 )) and all our data to 


f2{x) =Aq + ^ + ^ + ^ + ^. 

y/X X X'^'^ X^ 


(A.14) 


14a 


we 


The error err^rg is taken to be the maximum of the difference of dg ^ with these two other estimates. In fig. 
show the log]^o a function of Q for m/^ = 0.125, 0.3, 0.35, 0.5,1. It is clear that these errors are quite small. 

We have the largest error for mjg = 0.3 and Q = 0.5 which is explained by the fact that the gap is very small there 
as we are in the vicinity of a phase transition m- As mentioned above, it is well known that for smaller mass gaps, 
for a given bond dimension, the error on the ground-state MPS approximation will be larger. 

The continuum extrapolation of the electric field. 


^ = E (Lin)-Q)) (iV=|Z|), 

\nGZ[l,2Arl / ^ 


iez[i,2Ar] / g 

is found in a similar way. Now we use the values computed at x = 100, 200, 300,400 and perform a linear fit, 

gi{x) = Ag H— 

V X 


(A.15) 


(A.16) 


through the three largest x—values. The fact that we again have analytical behavior as a function of Xj^fx can be 
observed from fig. where we display the electric field as function of Ij^/x. It is also a consequence of the fact that 
{E{x)) = —daQ{x)/dQ and we already argued that cfq{x) is analytical as a function of x. To make our estim ate more 
robust against the choice of the interval and the fitting function we compute estimates by a linear fit (A.16) through 
all the points {x = 100, 200, 300,400) and a quadratic fit. 


g2{x) =Aq + ^ + 

^Jx X 


(A.17) 


through all the points. Again, the error err(^ ) is taken to be the maximum of the difference with these two estimates. 


The log]^o of orr^£;^ is displayed in fig. 14b 
transition at the critical mass (m/g) 


The errors are quite small but become larger again around the phase 
0.33 when going towards Q = 1/2. 
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FIG. 15: Q = 0.3. Continuum extrapolation of the electric field (E) for different values of m/g. 


At Q = 1/2 we do not display our error because this is a special case. For m/g < {m/g)c the CT symmetry is 
not broken and, thus, we shoul d have {E) = 0, and this for all values of x. Therefore a continuum extrapolation 
of {E) is useless, see fig. 16a To obtain a n err or bound we take the largest value in magnitude of {E{x)) for 

When m/g > {m/g\ we have two different vacua with opposite 


is useless, see fig. 1 
100,200,300,400. It is displayed in table III 


sign for the electric field. We will always take the negative sign which comes down to taking the vacuum in the limit 
Q ^ 1/2 for Q < 1/2. In this case it is possible to perform a polynomial extrapolation, see fig. 16b The results are 
given in Table [lll| If possible we compare with m- 
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0.4453 

-2.5507 

-5.5468 

-8.5429 

-11.5390 

-14.5351 


xIO'"^ 


0 


0.02 


{E[x)) [m/g = 0.3) 


O 


O 


o 

0.04 0.06 0.08 

l/v^ 

(a) 


O 


0.1 



FIG. 16: 


Q = 1/2. 


Continuum extrapolation of the electric field {E). (a): 


mj g = 0.3. 


(b): mjg = 0.35. 


m/g 

{E)/g 

{E)/g |60] 

0.125 

3 xl0“‘‘ 

- 

0.25 

2 xlO-^ 

- 

0.3 

0.0014 

0.0(3) 

0.35 

-0.313(2) 

- 

0.5 

-0.42041(3) 

-0.421(1) 

0.75 

-0.46145(2) 

- 

1 

-0.47692(2) 

-0.4769 (5) 

2 

-0.49364(3) 

- 

4 

-0.49834(3) 

- 


TABLE III: Electric field at Q = 1/2 for different value of m/^. 


3. Continuum extrapolation of the half-chain von Neumann entropy 


Using the Schmidt values see eq. (A.7), we can compute the half-chain von Neumann entropy Sq{x)^ 


Pmax 

Sq(x) = - E 

Q—Pmin CX.q — ^ 


for a particular value of x. As already mentioned in the main text, because the Schwinger model is equivalent to 
a noncritical boson theory mi, the half-chain von Neumann entropy should diverge as ( —1/6) log(l/y^) [87] when 
X Too. 

Let us check this for (5 = 0. Inspecting Sq{x) + ^ log(l/v^) as a function of x, we observe that it behaves linear 


m/g 

0.125 

0.25 

0.3 

0.35 

0.5 

0.75 

1 

Bo + 1/6 

1 

o 

T-1 

X 

T-1 

3 X 10“‘‘ 

5 X 10“^ 

5 X lO-^ 

8 X 10“^ 

1.3 X lO"'^ 

2 X 10“^ 


TABLE IV: The larg est val ue in magni tude of + 1/6 obtained from the fit (A. 18) through the largest five 
X—values and the fits (A. 18) and (A. 19) through all our data. According to |87f we should have 5o + 1/6 = 0. 
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FIG. 17: Q = 0. (a) m/^ = 0.125. Fit of the form ( — 1/6) log(l/v^) A-\- Cj-Jx through So{x). Inset: linear 

extrapolation of So{x) + (1/6) log(l/y^) based on the largest five x—values, x = 200, 300,400, 600, 800, to obtain the 
coefficients A and C. (b): same as (a) but now for m/g = 0.75. 


as a function of Xjy/x^ see inset figs. 17a and 17b Therefore, we should be able to fit So{x) to a function of the form 

Mx) = Ao + Bo\og(^^'^+Co^ (A.18) 

and find Bq = 1/6. Specifically, we fit our data corresponding to the largest five x—values, x = 200, 300,400, 600, 800, 
against fi to obtain a first estimate for Bq. To have some robustness against the choice of fitting interval and the 
fitting function, we also include our result for x = 100 and fitted all our data against /i and against 


f2{x) — Aq -\r Bq log ( ^^X 


(A.19) 


This gave us two other estimates for Bq. In table IV we give the results for Bq + 1/6. The value that is shown is 


the largest value for 5o + 1/6 (in magnitude) from the three fits, i.e. the largest error on the predicted result of 
l87l . As one observes these errors are at most 2 x 10“^ and for s mall values oi mjg only of order 10“^ which is a 


nice cross-check on our results. In the insets of figs. 


17a 


and 


17b 


we show a linear fit of the form f[x) = A + 
through Sq{x) + (l/6)log(x). Here we estimate A and C by taking into account the largest five x—values. In the 
main figures we also show the fit ( — 1/6) log(x) + f{x) through Sq{x). As expected, given the results in table IV, 
this fit matches our data very well. In fig. [^we show results for m/g = 0.125 and m/g = 0.75. A similar plot for 
m/g = 0.25 is shown in fig. 2a (main text). 


Because the coefficient of the logarithmic divergence of the von Neumann entropy is universal, the renormal¬ 
ized entropy AS'q = Sq — Sq should be UV-finite. In fig. 18 we plot ASq{x) as a function of l/^/x and observe that 
this scales linearly in l/^/x to the continuum limit. A continuum result for different values of Q and m/g is obtained 
in exactly the same way as for the electric field. The results are shown in the main text, see fig. [21^ The errors 
originating of the choice of fitting interval and fitting function are relatively small. 

For the electric field and the string tension we have that our results at x = 100, or equivalently ga = l/^/x = 0.1, 
only differ from the continuum result by at most 10 percent, see fig. and fig. In contrast, for the entropy 
this is not the case at all, see fig. the result at x = 100 and the continuum result differ by a factor of order one 
and sometimes also have a different sign. The main lesson is that, contrary to other quantities like the electric field 
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FIG. 18: Q = 0.3. Continuum extrapolation of the renormalized half-chain von Neumann entropy for different 
values of mjg. 


and the string tension, we should be careful when extrapolating results at finite x of the renormalized entropy to the 
continuum limit. In particular, for the nonuniform case, see main text section IV, one should always check how the 
results scale for different values of x. 
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FIG. 19: Some diagrams for the effective action from (B.l). On the first line we have the tree-level and the 
next-to-leading order jvr? contribution to Co (B.3). Evaluation of the first diagram on the second line would give 
a jm^ correction to Co, while the other diagram would give the leading g^jm^ contribution to C\. 


Appendix B: Perturbative calculation of gq 


To compute gq in the weak-coupling expansion we start from the Lagrangian (2.1) and include a current = ge^^d^Q] 
with Q constant everywhere in the bulk, and Q ^ 0 only at the boundaries at infinity (see my 


£ = ^ + gA^) 

= yp + gA^) - m) ^ - ^F^^F^" - ^F^.F'^" , (B.l) 


where on the last line we perform a partial integration and = e^^gQ. 

The effective action, obtained by integrating out both the fermion and the gauge fields in the path integral, will 
then have the general form: 

Seff= [d\Ceff= [d^x Co(-)F^.F^" + Ci(^)tS^^T!T + ... ^ (B.2) 

J J Tn TTh TYX/ 


where we can exclude derivative terms since F^j, is constant. At next-to-leading order we find for the first coefficient 
Co: 


C - b g' 

^ 4 247rm^ 


(B.3) 


The zero-order term here is the tree-level result while the g^ jw? term follows from the one-loop Feynman diagram 
on the first line of fig. 19, which can be calculated with standard techniques (see e.g. |lQ8j b Furthermore one can see 
that all other nonzero diagrams will lead to contributions to the coefficients Ci that are at least order g^jm^. Finally, 
we can then identify Seff = jdPxaq^ leading to the result (3.4). 


Appendix C: Details on the implementation of the DMRG optimization for the nonuniform case 


Consider the Schwinger Hamiltonian, 


H = 


^ (j2{L{n) + «(n)]2 + ^(-l)”(a,(n) + (-1)”) + x ^(cT+(n)F^(C^-(„ + i) + h.c.)) , (C.l) 

Vnez ^ nez nez y 
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in a nonuniform background field a{n) and assume we already computed a MPS approximation |^( 74 ( 1 ), ^4(2))) of 
the form (3.1) with virtual dimensions D{1) and D(2) for the zero-background Hamitonian {a{n) = 0) [57], i.e. 


|5'(A(l), A(2))) = ( n j Vr\k), 

hL \n£Z J 


with 


(C. 2 ) 


(C.3) 


Hn — ('^n 5 Pn)’) ^ { f 5 f I? Pn ^ ^[Pmin “1“ 1 ruod 2), Pmax “1“ 1 ruod 2)], | K-) — | {/^n}nGz) 5 

is a ground state of (0 with a{n) = 0. Note that the tensors A{n) take the form 

[^s,p('^)] ('g + i — = 1 • • • Dq(n mod 2), /3 = 1 . . . Drijl “h 1 mod 2),(C.4) 


qez 

G{n) 


Pmini'^ mod 2)^pmax{'^ mod 2)]; p^r G Z[p^^^(n + 1 mod 2),p^acc(^ + l mod 2)] in order to enforce Gauss’ law, 
4^(74(1), ^4(2))) = 0 with 


G(n) = L(n) - L(n - 1 ) - = q, 


(C.5) 


to the state (see eq. ( |2.9| )). We now consider a constant background electric field a{n) which has compact support: 
a{n) = a G M for n G N[0, /c], a{n) = 0 for n ^ N [Q, k]. The MPS trial state as ansatz for the ground state of this 
Hamiltonian that we will consider is, see eq. (4.1), 


^tr-I 




(C. 6 ) 


\n<rL 


in>rR 


where vl 0 < k vr and A^{n) = A^{n mod 2 ) corresponds to the MPS approximation ( |C. 2 [ ) of the ground state 
of the zero-backgro und H amiltonian. We take vr and vr odd. To enforce Gauss’ law, G{n) |4>(S)) = 0, the B{n) 
must take the form ( |G.4 ): 


[^s,p(^)]('g^CKq),(r,/ 3 r) ~ \Ps,p{l^)]^^^^^^q-\-{sn-\-i-l)^)/2,r^r,p^ aq — 1 . . . Dq{n)^ f3 — 1 . . . Dr{n + 1), 


(C.7) 


q e 'Z[Pmin{n),Pmax{n)]-, p,r E Z[pmin{n + '^),Pmax{n + 1)] where Dq(ri) = Dq{l), Dq{rR) = Dg{l), PminImaxirL) = 
Pmin/maxirR) = Prain/maxW- The formal virtual dimensions of this MPS are D{n) = YfqZltuA) Later in 

this appendix we will come back to the issue of which values to take for Dq{n) and Pmin/max{'^)' 


To obtain the best approximation within this class of states of the ground state of the Hamiltonian (G.l) we 
have to minimize 


H{b,b) = 


mB[b])\HmB[b])) 
($( b [ 5 ]) mB[b])) 


(C. 8 ) 


with respect to b{rL), • • •, b{rR — 1). This is a perfect problem to tackle with the DMRG [89]. We briefly sketch how 
this works in our case. 

The DMRG first minimizes H (6, b) with respect to b{rL) while keeping b{rL +1),..., b{rR — 1) fixed, then minimizes 
i7(5, b) with respect to 6(rL + 1) while keeping b{rL),b{rL + 2),..., b{rR — 1) fixed and so on until b{rR — 1). After 
this sweep, it will sweep back: minimizing 77(5, b) with respect to b{rR — 1) while keeping b{rL), • • •, b{rR — 2) fixed, 
then minimizing 77(5, b) with respect to b{rR — 2) while keeping b{rL), • • •, b{rR — 3), b{rR — 1) fixed and so on until 
b{rL)‘ The algorithm keeps sweeping until convergence of the quantity H{b^b) is reached. 

Let us now discuss how to minimise H{b^b) with respect to b{m) {vr ^ tu < vr — 1). It is convenient to use the 
gauge freedom m of the matrices: 

24,(1) ^ U{1)A^{1)U{2)-\A{2) ^ U{2)A^{2)U{l)-\B^{n) ^ V{n)B^{n)V{n+l)-\V{rL) = U{l),V{rR) = U{1), 

(C.9) 

to bring (|G.6|) in the following form 


' m—1 


MB)) = 534 I n n R I'") (^.10) 

Kn<rL / \n=rL / \n=m+l / \n>rR 
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where L{n), are in the left-canonical form: L,^{n) = 1, = 1, 

and Rn, B^^\n) are in the right-canonical form: ^ (Z]/^ = 

Z]s=-i 1 ^ ~ (<^ 5 ^))- III this case, the norm of the state is N[B{m)] = >/(^(5) |^(S)) = 

VE.tr[(B.( m)yB^{m)] which can be put to one by rescaling B{m): B{m) —> B{m)/N[B{m)]. Note that L{n) 
and R{n) depend only on the parity of n. Furthermore, we can use the remaining gauge freedom in the matrices to 
find positive-definite diagonal matrices A(n) with tr[A(n)] = 1, such that 


^ i:,(l)A(ri)(L.{l))t = A(ri - 1 ), ^ i.{2)A(ri - l)(i.(2))t = A(rz,), 

K K 

5;;Bl«(n)A(n)(B<«(n))* = A(» - l)(ri < „ <m-l) (C.lla) 


and 


5](i?(«)(n))tA(n - l)Bi^\n) = A(n)(m + 1 < n < r^), 5](i?,(l))tA(rfi - l)i?,(l) = Aim), 

K, K 

^(i?,(2))tA(rfl)i?,(2) = A(rfl + 1). (C.llb) 


Because the tensors R^L and B take the form (C.4) and (C.7), the diagonal elements of A^ can be labeled by the 
eigenvalues q of L{n)\ 

[-^{'^)][{qaq),{r/3r)] ~ ^Oiq,^r^q,r^q,ar{'^) 

Pmacc 

^q’) ^q 1 . . . Dq{Tl 1), 0 ^ ^g,Dq(n+l) {r) — '^g,Dq(n+l) —1 (^) ^ ^ ^ I5 ^ ^ ^ ^ ^q^aq I? 

(C.12) 

q^r e + l),Pmacc(^ + !)]• These diagonal elements of A(n) are the Schmidt values associated with the 

bipartition {Ai{n) = Z[— oo,n], ^ 2 (^) = + l,+oo]} of t he lat tice. More specifically, we have that the Schmidt 

decomposition with respect to this bipartition reads, see eq. (2. 


Pmax (n+l) Dq(n+1) 

mB))= Y1 E 

q=Pminin-\-l) Oiq = l 


(C.13) 


where (resp. )) are orthonormal unit vectors in the tensor product of the local Hilbert spaces in the 

region Ai{n) (resp. A 2 {ti)). At the boundaries (n < vl and n > vr) the Schmidt values depend only on the parity 
of n, more specifically: 


Xq,aq{'^ri - 1) = Xq,aq{rL) (Vn : 2n - 1 < tl), Xq,aq{‘^n) = Xq,aq{rL - l)(Vn : 2n < tl - 1,) 
Xq^aq (2^/ I) ~ Xq^(Xq {tr) (\/tI 2 tI 1 ^ , Xq^Q,^ (2'^) = ^q,0!.q ij'R T 1) (VtI : 2 tI < Tr -|- 1). 


(C.14a) 
(C.14b) 

These Schmidt values correspond to those of the gr ound state ( |C.2[ ) of the zero-background Hamiltonian and were 
already computed before (see appendix [A| subsection A1 for the details). 

With the MPS in the form ( |C.1Q[ ) Hih^h) is minimized with respect to bm by finding the smallest eigenvalue Eq 
and the corresponding eigenvector of the matrix H{m) with components 


[H[m)] (si,pi,ai,/3i);(s2,P2,CK2,/32) 


d 


d 


(^)]<ai/3i d[b S2,P2 i'rn)]a2P2 


{^im)\HmB[b])), 


(C.15) 


where € {-1,1},pfe € Z\pmin{m+l),Pmaxim+l)], Uk = 1... +(_i)fc)/ 2 (m),/3fc = 1... (m + 1) (k = 1,2). 

Because we are only interested in the smallest eigenvalue So and its eigenvector we can use the Lanczos iteration m- 
For thi s we only need the action of H{m) on h{m). Exploiting the gauge-invariant structure of the tensors R^L and 
H, see (C.4) and (C.7), the computation time of every sweep scales as 


Pmaxijl^ 

^1 E E ^D,in)f 

, n=rL-l q=pmin{n) 


o [ {vr - tl + 2 ) max(p max (tT/) Pmini'^)) 


m8iK{Dq{n)) 

n,q 
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We conclude this appendix by discussing how to fix Pminlmaxi'^) arid DgM - As choosing finite values for these 
quantities means an effective truncation in the Schmidt decomposition ( |C 13[ ) we need to look at the weight of the 
Schmidt values Xq^aqi'^) over the sectors q corresponding to the eigenvalues of L{n) for any n with < n < vr. 
Assuming that the ground-state approximations for n < vr and n > vr are accurate (see appendix [A| subsection [AT] ) , 
we don’t have to care about the Schmidt values (C.14) at the boundary. 

In practice we start with a certain distribution of values for each n, anticipating that the dominant eigenvalue 
sector of L{n) would shift from q = 0 at large n to g ^ Q at the center. After a first full DMRG-optimization, the 
initial Dq values are updated: increased in case that the minimal retained Schmidt value in the particular eigenvalue 
sector is larger than Xmin = decreased in case that the minimal retained Schmidt value is smaller. This is 

repeated a few times until all retained minimal Schmidt values are smaller or of the same order as Xmin- As for the 


choice of tl and tr^ we verify a posteriori that the inhomogeneous interval of the MPS (C.IO) is taken to be large 
enough, by verifying the convergence of local observables at large distances to their value for the homogeneous ground 
state. 

Let us give a specific example. In fig. and fig. we show some details of the simulation of the ground state 
for m/g = 0.25, Q = 5,x = 100, Lg = 10.1. In our setup with lattice spacing l|g‘^^/x = 0.1/g this corresponds 
to a distance of 101 sites between the external antiquark with charge —gQ and the external quark with charge 
gQ. Specifically, we put the antiquark at site 151 and the quark at site 252. And we reserve 150 sites on the left 
of the antiquark and 150 sites on the right for the nonuniform part of our MPS ansatz. In total we thus have 
151 + 101 150 = 402 tensors Bn that need to be optimized. By looking at the 10-base logarithm of the expectation 


value of some local quantities with respect to the Schwinger vacuum, see fig. |20a[ we observe that we take the range 
of the nonu niform part large enough: the errors by taking a finite range for the nonuniform part are of order 10“^. 
In fig. 20b we show the distribution of the minimum charge Pmin{'^) and maximum charge Pmax{'^) we used. For 


Q < Pniini^V and q > Pmaxi"^) we thus put Dq{n) = 0. The Pmini"^) and Pmaxi"^) we took at the boundaries, i.e. n > 1 
and n < 402 correspond to the Pmin and Pmax of the Schwinger vacuum, i.e. the vacuum without external charges, 
that we simulated in [57] . Between the boundaries and the external charges we anticipate the increasing electric field 
and raised Pmax{'^) to 4 + Q = 9 anticipating the dominant eigenvalue sector po ^ Q at the center. 

In figs. |20c| and |20d| we plot the distribution of the Schmidt values among the eigenvalues sector q of L{n) at the 
sites n = 150 (c) and n = 200 (d). As we explained above, we adapted the bond dimensio ns su ch t hat fo r each site 
n and at each eigenvalue sector p of L{n): miuc (n) < 10“^^. Comparing with figs. 


12a 


and 


12b 


we observe 


that the dominant eigenvalue sector is shifted to g = 2 for n = 150 and to g = 5 for n = 200. One can also see that 
our Pmini'^) ^ud Pmax{'^) are not entirely optimal: for certain charge sectors the largest Schmidt-val ue is still well 
below 10“^^, and these sectors could have been discarded altogether. As we can see by looking at fig. |21a| the most 
dominant eigenvalue sector of L{n), i.e. the eigenvalue sector g with the largest value for Ag,aq(^) shifts from 

g = 0 to g = 5 as we go from the left boundary to the middle and then decreases to g = 0 as we go to the right 
boundary. 

We also show the maximum bond dimension vaayiqDq(n) in fig. 
region where the electric background field is applied. 


21b The largest bond dimension is required in the 


Appendix D: The string state and broken-string state for Q = 1 in the weak-coupling limit 


As we mention in the introduction, a major difference between QCD and QED 2 is that the latter theory is already 
confining at the perturbative level, as the Coulomb potential is linear in 1+1 dimensions. This also allows us to 
understand the transition from the confining state to the broken-string two-meson state in the weak-coupling limit 
m/g ^ 1. In this nonrelativistic limit one can obtain the ground state by diagonalizing the Hamiltonian in subspaces 
of the different (fermion) particle sectors. The zero-particle sector simply consists of the Fock vacuum of the free 
Dirac field and corresponds to the confining string state with an energy 

£string= 9^L/2, (D.l) 


for probe charge Q = 1 and separation length L. The broken-string state will correspond to the ground state in the 
subspace of all states containing one (light) quark-antiquark pair. For this state the light antiquark will bind to the 
external probe quark and vice versa. We can make this more quantitative, by considering the effective Hamiltonian 
in the nonrelativistic limit for this particle sector: 


iJgg = 2 to - ^ ^ + y + L/2\ + —\xb- L/2\ + —\xa-xb\ - —\xa- L/2\ -—\xb + L/2\ + —L . (D.2) 
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(c) (d) 

FIG. 20: m/g = 0.25 ,x = 100, Q = b^Lg = 10.1. The stars represent the external charges. Between them the electric 
background field —Q = —5 is applied, (a): 10-base logarithm of the expectation values of some local quantities with 
the Schwinger vacuum-value subtracted. At the boundaries one observes that they are sufficiently small indicating 
that we took the nonuniform range wide enough, (b): maximum and minimum eigenvalues Pmax (^) and Pmin (^) of 
L{n — 1) we took into account in our numerical scheme on every site, (c): Distribution of the 10—base logarithm of 
the Schmidt values Xq^aq{i^) among the eigenvalue sectors q of L{n) for n = 150. (d): Distribution of the of the 
10—base logarithm of the Schmidt values Xq^aqi'^) among the eigenvalue sectors q of L{n) for n = 200. 
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(a) 

FIG. 21 : m/g = 0.25, x = 100, Q = 5, Lg = 10.1. The stars represent the external charges. Between them the electric 
background field —Q = —5 is applied.(a): Dominant eigenvalue sector of L{n), i.e. eigenvalue q of L{n) with largest 

Dq{n) is taken such that smallest Schmidt value is around 10 “^®. (b): maxqDq{n) : Largest bond 
dimension among the eigenvalue sectors of L{n) at every site n. 


Here xa and xb are the coordinates for the light antiquark and quark, and we put the probe quark at a; = —i /2 and 
the probe antiquark at a; = i/ 2 . Anticipating binding of the light fermions to the probe charges for large i, we can 
assume xa < xb, xa < i /2 and xb > —i /2 leading to a cancellation of the last four potential terms Hgq Ha + Hb 
with 

HA = m-^ + ^\xA + L/2\, 

HB = m-^ + ^\xB-L/2\. (D.3) 


A ground-state solution will therefore be of the form "^{xa^xb) = 4^a{xa)(I^b{xb) where now (I)a{xa) and (I)b{xb) 
are both ground states of the nonrelativistic one-particle problem for a linear potential. All eigenstates for this 
nonrelativistic Hamiltonian Ha (and similar for Hb) can be written in terms of the so-called Airy function Ai [96] : 

(Pa\^a) = MAi (^{g^my/^\xA + i/2| - , (D.4) 


where Af is the normalization factor and Sn is the (kinetic) eigenenergy of the eigenstate. These energies follow from the 
continuity requirement on <pA and ^(4 at xa = — A/ 2 , leading to either even or odd (pA under xa-\-L/2 ^ — (xa + A/ 2 ). 
The ground-state wave-function is even and the ground-state energy £0 is related to the first zero of the first derivative 

4/3 

of the Airy function, Ai'{xi) = 0 , at xi ^ —1.0188: Sq = —So in the nonrelativistic approximation we find: 


£2meson = 2 m + 1.0188 


g4/3 


m 


1/3 


(D.5) 


Notice that relativistic corrections to this approximation will necessarily involve quantum field contributions from 
other particle sectors. The relativistic one-particle Dirac equation has no bound state solutions for a linear (vector) 
potential nnuni]. 

Finally, in the nonrelativistic approximation we can then understand the transition from the string state to the 
broken-string state as a level crossing at the critical length L, where £ string = ^2meson- 











